library(plyr)
library(igraph)
library(RgoogleMaps)
library(RColorBrewer)
library(fmsb)
library(Hmisc)
library(MASS)
library(stargazer)
library(tidyverse)
library(multcomp)
library(WDI)
library(countrycode)
library(modelsummary) # https://modelsummary.com/articles/modelsummary.html#new-models-and-custom-statistics


setwd("~/Dropbox/China Huawei Paper/Production Materials/Replication/")
master.data <- read.csv("dataFile_China_Huawei.csv"); head(master.data)


# 250,000
dat <- master.data
dat$treated = 0
threshold=250000

dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(china.huawei > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

dat$ptALG = 1*(dat$cabb=="ALG")*dat$treated
dat$ptANG = 1*(dat$cabb=="ANG")*dat$treated
dat$ptBLR = 1*(dat$cabb=="BLR")*dat$treated
dat$ptBUI = 1*(dat$cabb=="BUI")*dat$treated
dat$ptCAM = 1*(dat$cabb=="CAM")*dat$treated
dat$ptCAO = 1*(dat$cabb=="CAO")*dat$treated
dat$ptCHA = 1*(dat$cabb=="CHA")*dat$treated
dat$ptDRC = 1*(dat$cabb=="DRC")*dat$treated
dat$ptDRV = 1*(dat$cabb=="DRV")*dat$treated
dat$ptETH = 1*(dat$cabb=="ETH")*dat$treated
dat$ptGAM = 1*(dat$cabb=="GAM")*dat$treated
dat$ptGUI = 1*(dat$cabb=="GUI")*dat$treated
dat$ptMYA = 1*(dat$cabb=="MYA")*dat$treated
dat$ptPAK = 1*(dat$cabb=="PAK")*dat$treated
dat$ptRWA = 1*(dat$cabb=="RWA")*dat$treated
dat$ptTAJ = 1*(dat$cabb=="TAJ")*dat$treated
dat$ptTUN = 1*(dat$cabb=="TUN")*dat$treated
dat$ptUGA = 1*(dat$cabb=="UGA")*dat$treated
dat$ptUKR = 1*(dat$cabb=="UKR")*dat$treated
dat$ptUZB = 1*(dat$cabb=="UZB")*dat$treated

dat = dat[dat$v2x_polyarchy<=0.42,]

m1 <- lm(v2smgovfilprc ~ ptALG + ptANG + ptBLR + ptBUI + ptCAM + ptCAO + ptCHA + ptDRC + ptDRV + ptETH + ptGAM + ptGUI + ptMYA + ptPAK + ptRWA + ptTAJ + ptTUN + ptUGA + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m1)

m2 <- lm(v2smgovshut ~ ptALG + ptANG + ptBLR + ptBUI + ptCAM + ptCAO + ptCHA + ptDRC + ptDRV + ptETH + ptGAM + ptGUI + ptMYA + ptPAK + ptRWA + ptTAJ + ptTUN + ptUGA + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m2)

m3 <- lm(v2smgovsmmon ~ ptALG + ptANG + ptBLR + ptBUI + ptCAM + ptCAO + ptCHA + ptDRC + ptDRV + ptETH + ptGAM + ptGUI + ptMYA + ptPAK + ptRWA + ptTAJ + ptTUN + ptUGA + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m3)

m4 <- lm(v2smarrest ~ ptALG + ptANG + ptBLR + ptBUI + ptCAM + ptCAO + ptCHA + ptDRC + ptDRV + ptETH + ptGAM + ptGUI + ptMYA + ptPAK + ptRWA + ptTAJ + ptTUN + ptUGA + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m4)

n.treated <- length(sort(unique(dat$cabb[dat$treated==1]))) - 2
cont <- matrix(c(0, sum(dat$ptALG,na.rm=T), sum(dat$ptANG,na.rm=T), sum(dat$ptBLR,na.rm=T), sum(dat$ptBUI,na.rm=T), sum(dat$ptCAM,na.rm=T), sum(dat$ptCAO,na.rm=T), sum(dat$ptCHA,na.rm=T), sum(dat$ptDRC,na.rm=T), sum(dat$ptDRV,na.rm=T), sum(dat$ptETH,na.rm=T), sum(dat$ptGAM,na.rm=T), sum(dat$ptGUI,na.rm=T), sum(dat$ptMYA,na.rm=T), sum(dat$ptPAK,na.rm=T), sum(dat$ptRWA,na.rm=T), sum(dat$ptTAJ,na.rm=T), sum(dat$ptTUN,na.rm=T), sum(dat$ptUGA,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptUZB,na.rm=T), rep(0, (length(coef(m4))-n.treated-1) )),1)

summary(glht(m1,cont))
summary(glht(m2,cont))
summary(glht(m3,cont))
summary(glht(m4,cont))

total.treatment.periods <- sum(sum(dat$ptALG,na.rm=T), sum(dat$ptANG,na.rm=T), sum(dat$ptBLR,na.rm=T), sum(dat$ptBUI,na.rm=T), sum(dat$ptCAM,na.rm=T), sum(dat$ptCAO,na.rm=T), sum(dat$ptCHA,na.rm=T), sum(dat$ptDRC,na.rm=T), sum(dat$ptDRV,na.rm=T), sum(dat$ptETH,na.rm=T), sum(dat$ptGAM,na.rm=T), sum(dat$ptGUI,na.rm=T), sum(dat$ptMYA,na.rm=T), sum(dat$ptPAK,na.rm=T), sum(dat$ptRWA,na.rm=T), sum(dat$ptTAJ,na.rm=T), sum(dat$ptTUN,na.rm=T), sum(dat$ptUGA,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptUZB,na.rm=T))
atreat.m1 <- glht(m1,cont/total.treatment.periods); summary(atreat.m1)
atreat.m2 <- glht(m2,cont/total.treatment.periods); summary(atreat.m2)
atreat.m3 <- glht(m3,cont/total.treatment.periods); summary(atreat.m3)
atreat.m4 <- glht(m4,cont/total.treatment.periods); summary(atreat.m4)

m1.250coef = summary(atreat.m1)$test[3]$coefficients 
m1.250se = summary(atreat.m1)$test[4]$sigma 
m1.250p = summary(atreat.m1)$test[6]$pvalues 

m2.250coef = summary(atreat.m2)$test[3]$coefficients 
m2.250se = summary(atreat.m2)$test[4]$sigma 
m2.250p = summary(atreat.m2)$test[6]$pvalues 

m3.250coef = summary(atreat.m3)$test[3]$coefficients 
m3.250se = summary(atreat.m3)$test[4]$sigma 
m3.250p = summary(atreat.m3)$test[6]$pvalues 

m4.250coef = summary(atreat.m4)$test[3]$coefficients 
m4.250se = summary(atreat.m4)$test[4]$sigma 
m4.250p = summary(atreat.m4)$test[6]$pvalues 

# 500,000
dat <- master.data
dat$treated = 0
threshold=500000

dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(china.huawei > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

dat$ptALG = 1*(dat$cabb=="ALG")*dat$treated
dat$ptANG = 1*(dat$cabb=="ANG")*dat$treated
dat$ptBUI = 1*(dat$cabb=="BUI")*dat$treated
dat$ptCAO = 1*(dat$cabb=="CAO")*dat$treated
dat$ptCHA = 1*(dat$cabb=="CHA")*dat$treated
dat$ptDRC = 1*(dat$cabb=="DRC")*dat$treated
dat$ptDRV = 1*(dat$cabb=="DRV")*dat$treated
dat$ptETH = 1*(dat$cabb=="ETH")*dat$treated
dat$ptGAM = 1*(dat$cabb=="GAM")*dat$treated
dat$ptGUI = 1*(dat$cabb=="GUI")*dat$treated
dat$ptMYA = 1*(dat$cabb=="MYA")*dat$treated
dat$ptPAK = 1*(dat$cabb=="PAK")*dat$treated
dat$ptRWA = 1*(dat$cabb=="RWA")*dat$treated
dat$ptTAJ = 1*(dat$cabb=="TAJ")*dat$treated
dat$ptTUN = 1*(dat$cabb=="TUN")*dat$treated
dat$ptUGA = 1*(dat$cabb=="UGA")*dat$treated
dat$ptUKR = 1*(dat$cabb=="UKR")*dat$treated
dat$ptUZB = 1*(dat$cabb=="UZB")*dat$treated

dat = dat[dat$v2x_polyarchy<=0.42,]

m1 <- lm(v2smgovfilprc ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptDRV + ptETH + ptGAM + ptGUI + ptMYA + ptPAK + ptRWA + ptTAJ + ptTUN + ptUGA + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m1)

m2 <- lm(v2smgovshut ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptDRV + ptETH + ptGAM + ptGUI + ptMYA + ptPAK + ptRWA + ptTAJ + ptTUN + ptUGA + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m2)

m3 <- lm(v2smgovsmmon ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptDRV + ptETH + ptGAM + ptGUI + ptMYA + ptPAK + ptRWA + ptTAJ + ptTUN + ptUGA + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m3)

m4 <- lm(v2smarrest ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptDRV + ptETH + ptGAM + ptGUI + ptMYA + ptPAK + ptRWA + ptTAJ + ptTUN + ptUGA + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m4)

n.treated <- length(sort(unique(dat$cabb[dat$treated==1]))) - 2
cont <- matrix(c(0, sum(dat$ptALG,na.rm=T), sum(dat$ptANG,na.rm=T), sum(dat$ptBUI,na.rm=T), sum(dat$ptCAO,na.rm=T), sum(dat$ptCHA,na.rm=T), sum(dat$ptDRC,na.rm=T), sum(dat$ptDRV,na.rm=T), sum(dat$ptETH,na.rm=T), sum(dat$ptGAM,na.rm=T), sum(dat$ptGUI,na.rm=T), sum(dat$ptMYA,na.rm=T), sum(dat$ptPAK,na.rm=T), sum(dat$ptRWA,na.rm=T), sum(dat$ptTAJ,na.rm=T), sum(dat$ptTUN,na.rm=T), sum(dat$ptUGA,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptUZB,na.rm=T), rep(0, (length(coef(m4))-n.treated-1) )),1)

summary(glht(m1,cont))
summary(glht(m2,cont))
summary(glht(m3,cont))
summary(glht(m4,cont))

total.treatment.periods <- sum(sum(dat$ptALG,na.rm=T), sum(dat$ptANG,na.rm=T), sum(dat$ptBUI,na.rm=T), sum(dat$ptCAO,na.rm=T), sum(dat$ptCHA,na.rm=T), sum(dat$ptDRC,na.rm=T), sum(dat$ptDRV,na.rm=T), sum(dat$ptETH,na.rm=T), sum(dat$ptGAM,na.rm=T), sum(dat$ptGUI,na.rm=T), sum(dat$ptMYA,na.rm=T), sum(dat$ptPAK,na.rm=T), sum(dat$ptRWA,na.rm=T), sum(dat$ptTAJ,na.rm=T), sum(dat$ptTUN,na.rm=T), sum(dat$ptUGA,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptUZB,na.rm=T))
atreat.m1 <- glht(m1,cont/total.treatment.periods); summary(atreat.m1)
atreat.m2 <- glht(m2,cont/total.treatment.periods); summary(atreat.m2)
atreat.m3 <- glht(m3,cont/total.treatment.periods); summary(atreat.m3)
atreat.m4 <- glht(m4,cont/total.treatment.periods); summary(atreat.m4)

m1.500coef = summary(atreat.m1)$test[3]$coefficients 
m1.500se = summary(atreat.m1)$test[4]$sigma 
m1.500p = summary(atreat.m1)$test[6]$pvalues 

m2.500coef = summary(atreat.m2)$test[3]$coefficients 
m2.500se = summary(atreat.m2)$test[4]$sigma 
m2.500p = summary(atreat.m2)$test[6]$pvalues 

m3.500coef = summary(atreat.m3)$test[3]$coefficients 
m3.500se = summary(atreat.m3)$test[4]$sigma 
m3.500p = summary(atreat.m3)$test[6]$pvalues 

m4.500coef = summary(atreat.m4)$test[3]$coefficients 
m4.500se = summary(atreat.m4)$test[4]$sigma 
m4.500p = summary(atreat.m4)$test[6]$pvalues 


# 1 million
dat <- master.data
dat$treated = 0
threshold=1000000

dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(china.huawei > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

dat$ptALG = 1*(dat$cabb=="ALG")*dat$treated
dat$ptANG = 1*(dat$cabb=="ANG")*dat$treated
dat$ptBUI = 1*(dat$cabb=="BUI")*dat$treated
dat$ptCAO = 1*(dat$cabb=="CAO")*dat$treated
dat$ptCHA = 1*(dat$cabb=="CHA")*dat$treated
dat$ptDRC = 1*(dat$cabb=="DRC")*dat$treated
dat$ptETH = 1*(dat$cabb=="ETH")*dat$treated
dat$ptGAM = 1*(dat$cabb=="GAM")*dat$treated
dat$ptGUI = 1*(dat$cabb=="GUI")*dat$treated
dat$ptPAK = 1*(dat$cabb=="PAK")*dat$treated
dat$ptTAJ = 1*(dat$cabb=="TAJ")*dat$treated
dat$ptTUN = 1*(dat$cabb=="TUN")*dat$treated
dat$ptUKR = 1*(dat$cabb=="UKR")*dat$treated
dat$ptUZB = 1*(dat$cabb=="UZB")*dat$treated

dat = dat[dat$v2x_polyarchy<=0.42,]

m1 <- lm(v2smgovfilprc ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptGAM + ptGUI + ptPAK + ptTAJ + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m1)

m2 <- lm(v2smgovshut ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptGAM + ptGUI + ptPAK + ptTAJ + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m2)

m3 <- lm(v2smgovsmmon ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptGAM + ptGUI + ptPAK + ptTAJ + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m3)

m4 <- lm(v2smarrest ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptGAM + ptGUI + ptPAK + ptTAJ + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m4)

n.treated <- length(sort(unique(dat$cabb[dat$treated==1]))) - 2
cont <- matrix(c(0, sum(dat$ptALG,na.rm=T), sum(dat$ptANG,na.rm=T), sum(dat$ptBUI,na.rm=T), sum(dat$ptCAO,na.rm=T), sum(dat$ptCHA,na.rm=T), sum(dat$ptDRC,na.rm=T), sum(dat$ptETH,na.rm=T), sum(dat$ptGAM,na.rm=T), sum(dat$ptGUI,na.rm=T), sum(dat$ptPAK,na.rm=T), sum(dat$ptTAJ,na.rm=T), sum(dat$ptTUN,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptUZB,na.rm=T), rep(0, (length(coef(m4))-n.treated-1) )),1)

summary(glht(m1,cont))
summary(glht(m2,cont))
summary(glht(m3,cont))
summary(glht(m4,cont))

total.treatment.periods <- sum(sum(dat$ptALG,na.rm=T), sum(dat$ptANG,na.rm=T), sum(dat$ptBUI,na.rm=T), sum(dat$ptCAO,na.rm=T), sum(dat$ptCDI,na.rm=T), sum(dat$ptCHA,na.rm=T), sum(dat$ptDRC,na.rm=T), sum(dat$ptETH,na.rm=T), sum(dat$ptGAM,na.rm=T), sum(dat$ptPAK,na.rm=T), sum(dat$ptTAJ,na.rm=T), sum(dat$TUN,na.rm=T), sum(dat$ptUZB,na.rm=T))
atreat.m1 <- glht(m1,cont/total.treatment.periods); summary(atreat.m1)
atreat.m2 <- glht(m2,cont/total.treatment.periods); summary(atreat.m2)
atreat.m3 <- glht(m3,cont/total.treatment.periods); summary(atreat.m3)
atreat.m4 <- glht(m4,cont/total.treatment.periods); summary(atreat.m4)

m1.1mcoef = summary(atreat.m1)$test[3]$coefficients 
m1.1mse = summary(atreat.m1)$test[4]$sigma 
m1.1mp = summary(atreat.m1)$test[6]$pvalues 

m2.1mcoef = summary(atreat.m2)$test[3]$coefficients 
m2.1mse = summary(atreat.m2)$test[4]$sigma 
m2.1mp = summary(atreat.m2)$test[6]$pvalues 

m3.1mcoef = summary(atreat.m3)$test[3]$coefficients 
m3.1mse = summary(atreat.m3)$test[4]$sigma 
m3.1mp = summary(atreat.m3)$test[6]$pvalues 

m4.1mcoef = summary(atreat.m4)$test[3]$coefficients 
m4.1mse = summary(atreat.m4)$test[4]$sigma 
m4.1mp = summary(atreat.m4)$test[6]$pvalues 


# 5 million
dat <- master.data
dat$treated = 0
threshold=5000000

dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(china.huawei > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

dat$ptALG = 1*(dat$cabb=="ALG")*dat$treated
dat$ptANG = 1*(dat$cabb=="ANG")*dat$treated
dat$ptBUI = 1*(dat$cabb=="BUI")*dat$treated
dat$ptCAO = 1*(dat$cabb=="CAO")*dat$treated
dat$ptCHA = 1*(dat$cabb=="CHA")*dat$treated
dat$ptDRC = 1*(dat$cabb=="DRC")*dat$treated
dat$ptETH = 1*(dat$cabb=="ETH")*dat$treated
dat$ptGAM = 1*(dat$cabb=="GAM")*dat$treated
dat$ptPAK = 1*(dat$cabb=="PAK")*dat$treated
dat$ptTAJ = 1*(dat$cabb=="TAJ")*dat$treated
dat$ptTUN = 1*(dat$cabb=="TUN")*dat$treated
dat$ptUKR = 1*(dat$cabb=="UKR")*dat$treated
dat$ptUZB = 1*(dat$cabb=="UZB")*dat$treated

dat = dat[dat$v2x_polyarchy<=0.42,]

m1 <- lm(v2smgovfilprc ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptGAM + ptPAK + ptTAJ + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m1)

m2 <- lm(v2smgovshut ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptGAM + ptPAK + ptTAJ + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc  + as.factor(cabb) + as.factor(year), data=dat) 
summary(m2)

m3 <- lm(v2smgovsmmon ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptGAM + ptPAK + ptTAJ + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc  + as.factor(cabb) + as.factor(year), data=dat) 
summary(m3)

m4 <- lm(v2smarrest ~ ptALG + ptANG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptGAM + ptPAK + ptTAJ + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m4)

n.treated <- length(sort(unique(dat$cabb[dat$treated==1]))) - 2
cont <- matrix(c(0, sum(dat$ptALG,na.rm=T), sum(dat$ptANG,na.rm=T), sum(dat$ptBUI,na.rm=T), sum(dat$ptCAO,na.rm=T), sum(dat$ptCHA,na.rm=T), sum(dat$ptDRC,na.rm=T), sum(dat$ptETH,na.rm=T), sum(dat$ptGAM,na.rm=T), sum(dat$ptPAK,na.rm=T), sum(dat$ptTAJ,na.rm=T), sum(dat$ptTUN,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptUZB,na.rm=T), rep(0, (length(coef(m4))-n.treated-1) )),1)

summary(glht(m1,cont))
summary(glht(m2,cont))
summary(glht(m3,cont))
summary(glht(m4,cont))

total.treatment.periods <- sum(sum(dat$ptALG,na.rm=T), sum(dat$ptANG,na.rm=T), sum(dat$ptBUI,na.rm=T), sum(dat$ptCAO,na.rm=T), sum(dat$ptCHA,na.rm=T), sum(dat$ptDRC,na.rm=T), sum(dat$ptETH,na.rm=T), sum(dat$ptGAM,na.rm=T), sum(dat$ptPAK,na.rm=T), sum(dat$ptTAJ,na.rm=T), sum(dat$ptTUN,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptUZB,na.rm=T))
atreat.m1 <- glht(m1,cont/total.treatment.periods); summary(atreat.m1)
atreat.m2 <- glht(m2,cont/total.treatment.periods); summary(atreat.m2)
atreat.m3 <- glht(m3,cont/total.treatment.periods); summary(atreat.m3)
atreat.m4 <- glht(m4,cont/total.treatment.periods); summary(atreat.m4)

m1.5mcoef = summary(atreat.m1)$test[3]$coefficients 
m1.5mse = summary(atreat.m1)$test[4]$sigma 
m1.5mp = summary(atreat.m1)$test[6]$pvalues 

m2.5mcoef = summary(atreat.m2)$test[3]$coefficients 
m2.5mse = summary(atreat.m2)$test[4]$sigma 
m2.5mp = summary(atreat.m2)$test[6]$pvalues 

m3.5mcoef = summary(atreat.m3)$test[3]$coefficients 
m3.5mse = summary(atreat.m3)$test[4]$sigma 
m3.5mp = summary(atreat.m3)$test[6]$pvalues 

m4.5mcoef = summary(atreat.m4)$test[3]$coefficients 
m4.5mse = summary(atreat.m4)$test[4]$sigma 
m4.5mp = summary(atreat.m4)$test[6]$pvalues 


# 10 million
dat <- master.data
dat$treated = 0
threshold=10000000

dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(china.huawei > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

dat$ptALG = 1*(dat$cabb=="ALG")*dat$treated
dat$ptBUI = 1*(dat$cabb=="BUI")*dat$treated
dat$ptCAO = 1*(dat$cabb=="CAO")*dat$treated
dat$ptCHA = 1*(dat$cabb=="CHA")*dat$treated
dat$ptDRC = 1*(dat$cabb=="DRC")*dat$treated
dat$ptETH = 1*(dat$cabb=="ETH")*dat$treated
dat$ptTUN = 1*(dat$cabb=="TUN")*dat$treated
dat$ptUKR = 1*(dat$cabb=="UKR")*dat$treated
dat$ptUZB = 1*(dat$cabb=="UZB")*dat$treated

dat = dat[dat$v2x_polyarchy<=0.42,]

m1 <- lm(v2smgovfilprc ~ ptALG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m1)

m2 <- lm(v2smgovshut ~ ptALG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc  + as.factor(cabb) + as.factor(year), data=dat) 
summary(m2)

m3 <- lm(v2smgovsmmon ~ ptALG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc  + as.factor(cabb) + as.factor(year), data=dat) 
summary(m3)

m4 <- lm(v2smarrest ~ ptALG + ptBUI + ptCAO + ptCHA + ptDRC + ptETH + ptTUN + ptUKR + ptUZB + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc + as.factor(cabb) + as.factor(year), data=dat) 
summary(m4)

n.treated <- length(sort(unique(dat$cabb[dat$treated==1]))) - 1
cont <- matrix(c(0, sum(dat$ptALG,na.rm=T), sum(dat$ptBUI,na.rm=T), sum(dat$ptCAO,na.rm=T), sum(dat$ptCHA,na.rm=T), sum(dat$ptDRC,na.rm=T), sum(dat$ptETH,na.rm=T), sum(dat$ptTUN,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptUZB,na.rm=T), rep(0, (length(coef(m4))-n.treated-1) )),1)

summary(glht(m1,cont))
summary(glht(m2,cont))
summary(glht(m3,cont))
summary(glht(m4,cont))

total.treatment.periods <- sum(sum(dat$ptALG,na.rm=T), sum(dat$ptBUI,na.rm=T), sum(dat$ptCAO,na.rm=T), sum(dat$ptCHA,na.rm=T), sum(dat$ptDRC,na.rm=T), sum(dat$ptETH,na.rm=T), sum(dat$ptTUN,na.rm=T), sum(dat$ptUKR,na.rm=T), sum(dat$ptUZB,na.rm=T))
atreat.m1 <- glht(m1,cont/total.treatment.periods); summary(atreat.m1)
atreat.m2 <- glht(m2,cont/total.treatment.periods); summary(atreat.m2)
atreat.m3 <- glht(m3,cont/total.treatment.periods); summary(atreat.m3)
atreat.m4 <- glht(m4,cont/total.treatment.periods); summary(atreat.m4)

m1.10mcoef = summary(atreat.m1)$test[3]$coefficients 
m1.10mse = summary(atreat.m1)$test[4]$sigma 
m1.10mp = summary(atreat.m1)$test[6]$pvalues 

m2.10mcoef = summary(atreat.m2)$test[3]$coefficients 
m2.10mse = summary(atreat.m2)$test[4]$sigma 
m2.10mp = summary(atreat.m2)$test[6]$pvalues 

m3.10mcoef = summary(atreat.m3)$test[3]$coefficients 
m3.10mse = summary(atreat.m3)$test[4]$sigma 
m3.10mp = summary(atreat.m3)$test[6]$pvalues 

m4.10mcoef = summary(atreat.m4)$test[3]$coefficients 
m4.10mse = summary(atreat.m4)$test[4]$sigma 
m4.10mp = summary(atreat.m4)$test[6]$pvalues 


# Export table

ti = data.frame(rbind(c(m1.250coef, m1.250se, m1.250p),c(m1.500coef, m1.500se, m1.500p),c(m1.1mcoef, m1.1mse, m1.1mp),c(m1.5mcoef, m1.5mse, m1.5mp),c(m1.10mcoef, m1.10mse, m1.10mp)))
ti$term = c("Transfer \u2265 $250,000","Transfer \u2265 $500,000","Transfer \u2265 $1,000,000","Transfer \u2265 $5,000,000","Transfer \u2265 $10,000,000")
colnames(ti) = c("estimate","std.error","p.value","term")
ti = ti[c("term","estimate","std.error","p.value")]
out1 = list(tidy=ti)
class(out1) = "modelsummary_list"

ti = data.frame(rbind(c(m2.250coef, m2.250se, m2.250p), c(m2.500coef, m2.500se, m2.500p), c(m2.1mcoef, m2.1mse, m2.1mp), c(m2.5mcoef, m2.5mse, m2.5mp),c(m2.10mcoef, m2.10mse, m2.10mp)))
ti$term = c("Transfer \u2265 $250,000","Transfer \u2265 $500,000","Transfer \u2265 $1,000,000","Transfer \u2265 $5,000,000","Transfer \u2265 $10,000,000")
colnames(ti) = c("estimate","std.error","p.value","term")
ti = ti[c("term","estimate","std.error","p.value")]
out2 = list(tidy=ti)
class(out2) = "modelsummary_list"


ti = data.frame(rbind(c(m3.250coef, m3.250se, m3.250p),c(m3.500coef, m3.500se, m3.500p),c(m3.1mcoef, m3.1mse, m3.1mp),c(m3.5mcoef, m3.5mse, m3.5mp),c(m3.10mcoef, m3.10mse, m3.10mp)))
ti$term = c("Transfer \u2265 $250,000","Transfer \u2265 $500,000","Transfer \u2265 $1,000,000","Transfer \u2265 $5,000,000","Transfer \u2265 $10,000,000")
colnames(ti) = c("estimate","std.error","p.value","term")
ti = ti[c("term","estimate","std.error","p.value")]
out3 = list(tidy=ti)
class(out3) = "modelsummary_list"

ti = data.frame(rbind(c(m4.250coef, m4.250se, m4.250p),c(m4.500coef, m4.500se, m4.500p),c(m4.1mcoef, m4.1mse, m4.1mp),c(m4.5mcoef, m4.5mse, m4.5mp),c(m4.10mcoef, m4.10mse, m4.10mp)))
ti$term = c("Transfer \u2265 $250,000","Transfer \u2265 $500,000","Transfer \u2265 $1,000,000","Transfer \u2265 $5,000,000","Transfer \u2265 $10,000,000")
colnames(ti) = c("estimate","std.error","p.value","term")
ti = ti[c("term","estimate","std.error","p.value")]
out4 = list(tidy=ti)
class(out4) = "modelsummary_list"


myrows = data.frame(list(c("Country Fixed Effects",rep("Yes",4))
                         ,c("Year Fixed Effects",rep("Yes",4))
                         ,c("Control Variables",rep("Yes",4))
))

export = modelsummary(list("Internet Filtering" = out1
                           , "Internet Shutdowns" = out2
                           ,"Social Media Monitoring"=out3
                           ,"Arrests for Political Content"=out4
)
,add_rows=data.frame(t(myrows))
, stars=c('*' = .1, '**' = 0.05, '***' = .01)
, output="latex_tabular"
)
export = gsub("≥", "$\\\\ge$", export)
export

